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A  DYNAMIC  MODEL  OF  AN  AXISYMMETRIC,  TRANSVERSELY  ISOTROPIC, 
FLUID-LOADED,  FULLY  ELASTIC  CYLINDRICAL  SHELL 


1.  INTRODUCTION 

Fluid-loaded  shells  are  encountered  in  numerous  natural  and  man-made  applications. 
Examples  include  arteries,  inner  ear  tubes,  hydraulic  lines,  marine  pilings,  water  pipes,  and  shock 
absorbers.  Understanding  the  behavior  of  these  systems  is  important  so  that  their  performance 
can  be  analyzed  or,  in  the  case  of  mechanical  systems,  the  next  generation  can  be  better 
designed.  Extensive  modeling  of  these  systems  has  been  done  over  the  years,  and  a  large  volume 
of  research  articles  exists  in  the  area  of  cylindrical  and  spherical  shells.  From  a  complexity 
standpoint,  membrane  models  of  shells  are  the  most  simplistic  and  have  been  derived  and 
analyzed  most  notably  by  Love.'  A  bending  stiffness  term  was  added  to  the  membrane  equations 
by  Donnell.  Rotary  inertia  and  shear  effects  were  added  to  the  cylindrical  shell  equations  by 
Mirsky  and  Herrmann.^  ”’  These  previous  models  are  based  on  membrane  and  flexural  wave 
theory  and  are  accurate  only  at  low  frequencies  and  low  wavenumbers.  The  fully  elastic, 
isotropic  cylindrical  shell  was  modeled  and  analyzed  by  Gazis.^’^  The  work  by  Gazis  was  a 
significant  extension  of  previous  theory  as  it  allowed  analytical  modeling  at  all  frequencies  and 
wavenumbers,  rather  than  just  a  small  subset  of  low-frequency  and  low- wavenumber  analysis. 
The  fiilly  elastic,  transversely  isotropic  cylindrical  shell  was  modeled  by  Laverty^  for  analysis  of 
wood  cylinders.  Fay^  added  fluid  loading  to  the  membrane  theories  for  a  solid  cylinder. 
Additionally,  Peloquin^  added  fluid  loading  to  various  flexural  wave  theories  for  hollow 
cylinders. 

This  report  develops  an  analytical  model  of  a  transversely  isotropic,  fluid-loaded,  fully 
elastic  axisymmetric  cylinder  that  is  in  contact  with  fluid  on  both  its  interior  and  exterior.  The 
model  begins  with  the  equations  of  motion  of  a  transversely  isotropic  body  in  cylindrical 
coordinates.  Using  the  radial  and  longitudinal  equations  of  motion,  two  free  wavenumbers  are 
calculated  corresponding  to  two  specific  waves  that  are  propagating  in  the  medium.  A  solution 
set  to  the  shell  displacement  field  is  formulated  that  contains  four  unknown  wave  propagation 
coefficients.  These  coefficients  are  inserted  into  the  stress  boundary  conditions  at  the  inner  and 


outer  surfaces  of  the  shell.  Also  included  in  these  boundary  conditions  are  the  pressure  loads  of 
the  inner  and  outer  fluid  fields  and  any  external  loads  that  may  be  acting  on  the  system.  This 
produces  four  algebraic  equations  with  four  unknown  wave  propagation  coefficients.  This  set  of 
equations  can  be  solved  to  obtain  an  analytical  solution  to  the  shell  displacements,  the  pressure 
of  the  inner  fluid,  and  the  pressure  of  the  outer  fluid.  The  model  is  verified  by  comparing  the 
results  with  two  previously  derived  models,  and  a  numerical  example  is  included  to  illustrate  the 
behavior  of  a  thick  shell  under  two  loading  conditions.  Additionally,  a  MATLAB  subroutine  is 
included  that  contains  a  vectorized  computation  that  outputs  interior  shell  pressure  produced 
from  external  forcing  functions. 
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2.  SYSTEM  MODEL 


The  system  equations  consist  of  three  separate  models:  the  cylindrical  shell  equations  of 
motion  in  the  radial  and  axial  direction,  the  inner  acoustic  field  wave  equation  of  pressure,  and 
the  outer  acoustic  field  wave  equation  of  pressure.  Once  the  general  solutions  to  these  equations 
of  motions  and  pressure  are  determined,  they  are  coupled  using  linear  momentum  and  inserted 
into  the  stress  fields  at  the  inner  and  outer  radii  of  the  shell.  This  produces  a  four-by-four  matrix 
that  contains  the  dynamics  of  the  system  multiplied  by  a  four-by-one  vector  that  contains  the 
unknown  wave  propagation  coefficients  and  is  equal  to  a  four-by-one  vector  containing  the 
applied  external  loads.  This  matrix  equation  can  be  solved  and  the  response  of  the  system  can  be 
calculated.  This  process  is  described  below.  A  schematic  of  the  system  illustrating  the 
coordinate  system  is  shown  in  figure  1 . 
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Figure  1.  Fluid-Loaded  Shell  with  Coordinate  System 
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The  equation  of  motion  of  a  fully  elastic,  isotropic  body  in  cylindrical  coordinates’®  in  the 
radial  direction  is 


d^u(x,r,t) 


d^u(x,rj)  1  du{x,r,t)  u{x,r,t) 
dr^  r  dr 
d^w{x,r,t) 


+  c 


44 


d^u{x,r,t) 


+  (c,3+C44) 


drdx 


(1) 


in  the  longitudinal  direction,  the  equation  is 


P- 


d^w{x,r,t) 
dt^ 


-  (^13  +C44) 


d^u{x,r,t)  ^  1  du{x,r,t) 


drdx 


dx 


+  C44 


d^w(x,r,t)  ^  1  d\v{x,r,t) 


dr^ 
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d  w{x,r,t) 
+  C33 - - - 

dx^ 


(2) 


In  equations  (1)  and  (2),  u{x,r,t)  is  the  displacement  in  the  radial  direction  (m),  w{x,r,t)  is  the 
displacement  in  the  longitudinal  direction  (m),  r  is  the  coordinate  of  the  radial  direction  (m),  x  is 
the  coordinate  of  the  longitudinal  direction  (m),  t  is  time  (s),  p  is  density  of  the  shell  (kg  m'^),  and 
Cij  are  stiffness  constants  that  contain  the  material  properties  (N  m‘^)  and  are  typically  complex 
quantities.  These  constants  are  determined  using  the  constitutive  equations  in  cylindrical 
coordinates  between  strain  and  stress  written  for  a  solid  that  is  transversely  isotropic  in  the  axial 
direction  with  respect  to  the  radial  and  circumferential  directions.  These  equations  are 
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Yxr  = 
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(7) 
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jcr  ’ 


and 


YxO  = 


(8) 


where  e\\  are  the  normal  strains  (dimensionless),  /ij  are  the  shear  strains  (dimensionless),  ca  are 
the  normal  stresses  (N  m'^),  rij  are  the  shear  stresses  (N  m'^),  Er  is  Young’s  modulus  in  the  radial 
direetion  (N  m'  ),  Ex  is  Young’s  modulus  in  the  axial  direction  (N  m' ),  Urx  is  Poisson’s  ratio  in 
the  longitudinal  direction  with  a  load  being  applied  in  the  radial  direction  (dimensionless),  and 
Uxr  is  Poisson’s  ratio  in  the  radial  direction  with  a  load  being  applied  in  the  longitudinal 
direction  (dimensionless).  Equations  (3)  through  (8)  are  inverted  so  that  the  stresses  are 
functions  of  the  strains,  which  in  matrix  form  is 


a  =  C^  , 


(9) 


where 


^  ”  [p^rr  ^90  ^xx  ^ xO  ^ xr  ^ rO  ]  ’ 


(10) 


and 


Z  £q0  YxO  Yxr  Yre\ 


(11) 


Using  equation  (9)  and  Betti’s  reciprocal  law  for  composite  materials,  written  as 


^rx  _  ^xr 


Er 


(12) 


the  stiffness  constants  in  equations  (1)  and  (2)  can  be  solved  for  in  terms  of  engineering 
constants.  They  are 
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(14) 


<^12 
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EfUxx 
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C33  = 


E,Uxr{\-Orx) 

^^rxi^-^^rx-^^rx^xr) 
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C44  -  Gxr  - 


2(l  +  f^;tr) 


(15) 

(16) 

(17) 


The  solution  to  equations  (1)  and  (2)  is  now  determined  for  free  wave  propagation  in  a 
medium  that  is  bounded  in  the  radial  direction,  unbounded  in  the  axial  direction,  and  harmonic  in 
time.  The  argument  is  made^  that  the  solution  to  the  transversely  isotropic  differential  equations 
has  to  have  the  same  form  as  the  solution  to  the  isotropic  differential  equations.  Thus,  the 
solution  in  the  radial  direction  is  written  as 


u(r,x,t)  =  C/(r)  exp(iAx)  exp(-i  6>r)  =  GBi(;r)exp(iAx)exp(-i(y/)  ,  (18) 

and  the  solution  in  the  longitudinal  direction  is  written  as 

w{r,x,t)  =  W{r) exp(ifcc) exp(-i o)t)  =  exp(ifcc) exp(- i oyt) ,  (19) 

where  G  and //are  unknown  wave  propagation  coefficients,  B]  denotes  an  ordinary  Bessel 
function  of  order  one,  Bq  denotes  an  ordinary  Bessel  function  of  order  zero,  /  is  the  free 
propagation  wavenumber  (rad  m  '),  k  is  the  wavenumber  with  respect  to  the  x-axis  (rad  m'*),  and 
i  is  .  It  is  noted  that  the  free  propagation  wavenumbers  are  typically  complex  quantities; 
thus,  the  Bessel  functions  contain  complex  arguments.  To  facilitate  this  type  of  analysis,  the 
Bessel  functions  will  all  be  ordinary  Bessel  functions  of  the  first  and  second  kind  with  complex 
arguments,  rather  than  switching  between  normal  and  modified  Bessel  functions  based  on  the 
sign  of  the  argument.  Substituting  equations  (18)  and  (19)  into  equations  (1)  and  (2)  yields  the 
two-by-two  system  of  algebraic  equations,  written  as 
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(20) 


2  2  2 
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+C44) 
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o' 
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H 

0 

The  determinant  of  the  two-by-two  matrix  in  equation  (20)  must  be  zero  if  a  solution  other  than 
the  trivial  solution  is  going  to  exist.  This  yields  a  quadratic  equation  with  respect  to  the 
propagation  wavenumber  that  is  written  as 

+  hy^  +  c  =  0  ,  (21) 

where 

a-CiiC44  ,  (22) 

b  =  (c\\C23  -C\2-2ci2C44)k^  -{C44  +Cii)pco^  ,  (23) 

and 

c  =  p^co^  -  (C33  +  C44  )poP'k^  +  033044^:^  .  (24) 


The  solution  to  equation  (21)  is 


y\,2 


-b±{b'^ -Aac)^'^ 
2a 


(25) 


Only  the  positive  values  from  equation  (25)  are  needed,  as  the  zero-order  and  first-order  Bessel 
functions  in  equations  (18)  and  (19)  are  even  functions;  thus,  negative  values  will  not  contribute 
to  a  linearly  independent  solution.  The  first  row  of  equation  (20)  yields 


^  _  pco^  -c^k^  ^ 

iA:/,  2(0,3+044) 


(26) 
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The  solution  is  now  written  as  Bessel  functions  of  the  first  kind  using  the  wavenumbers  y\ 
and  /2-  The  expressions  for  the  displacement  fields  are 

m(x,  r,  0  =  [G|  J I  (r,  /-)  +  G3 J ,  (r2  '■)]  exp(i Ax)  exp(-i  cot) ,  (2  7) 


H<x,  r,  /)  =  [//,  Jo  (/,  r)  +  7/3  Jo  (r^r)]  exp(i  Ax)  exp(-i  cot) 

=  [Gi^ihir\f)  +  Gi<^2hir2f)] exp(i Ax) exp(-i cot) . 


(28) 


Additionally,  because  the  domain  of  r  is  from  a{>  0)  to  b{<  co) ,  Bessel  functions  of  the  second 
kind  are  admissible  solutions,  and  the  expressions  for  the  displacement  fields  using  these 
functions  are 


M(x,r,/)  =  [G2Y,(y,r)  +  G4Y,(/2r)]exp(iAx)exp(-i(y/) , 

and 


(29) 


w{x,  r,  /)  =  [//j  Yo  (y,  r)  +  (/jr)]  exp(i  Ax)  exp(-i  cot) 

=  [<^2^1  Yo(/iO  +  <^4^2Yo(r2'')]exp(iAx)exp(-i(y/) . 


(30) 


The  problem  set  represented  by  equations  (27)  and  (28)  is  linearly  independent  from  the  solution 
set  given  by  equations  (29)  and  (30),  and  a  complete  solution  is  a  linear  combination  of  both 
equation  sets.  This  gives  the  total  solution  to  the  shell  displacements  as 

m(x, r,  0  =  [<7,  J ,  (y,r)  +  G2  Y,  (y,  r)  +  Gj J ,  (y^r)  +  G4 Y,  (y^r)]  exp(iAx)  exp(-icy 0  ,  (31) 

and 


H<x,  r,  /)  =  [G,  J  0  (/,  r)  +  G,  I,  Y^  (/,/-)  +  Gj  0  (^2 '')  +  <^4  ^2^0  (^2  exp(iAx)  exp(-i  cot),  (32) 

where  Gj,  G2,  G3,  and  G4  are  unknown  wave  propagation  coefficients.  The  insertion  of 
equations  (31)  and  (32)  into  equations  (1)  and  (2)  verifies  that  they  are  solutions  to  the  original 
differential  equations  of  motion. 
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The  unknown  wave  propagation  coefficients  are  determined  using  the  four  stress-boundary 
conditions  of  the  shell.  The  first  boundary  condition  is  a  force  balance  between  the  pressure  in 
the  interior  fluid  and  the  normal  radial  stress  in  the  shell  at  the  interface  where  r=  a.  This 
equation  is  written  as 


du(x,a,t)  u{x,a,t)  dw{x,a,t) 

+  Ci2 - +  ci3 - - =  -pi{x,a,t), 


dr 


a 


dx 


(33) 


where  p^{x,a,t)  is  the  pressure  of  the  interior  fluid  (N  m'^)  at  r  =  a,  which  satisfies  the  wave 
equation  in  cylindrical  coordinates;  i.e., 


d^ Pi{x,r,t)  ^  1  dpi{x,r,t)  ^  d^ p;(x,r,t) _ 1  d'^ pi(x,r,t)  ^  ^ 

dr^  r  Sr  cf  dt^ 

where  C/  is  the  acoustic  (or  compressional)  wavespeed  of  the  interior  fluid  (m  s'').  Using  the 
infinite  length  of  the  cylinder  in  the  jc-direction  and  the  constraint  that  the  pressure  field  has  to  be 
finite  at  r  =  0,  the  temporal  harmonic  solution  to  equation  (34)  is 

Pi(x,r,t)  =  /^(r)exp(iAx)exp(-i<y/)  =  A/Jo(y,r)exp(ifcc)exp(-i<w/) ,  (35) 

where 


Vi  = 


\ 

CO 


H  1/2 


-k^ 


(36) 


In  general,  y,  can  be  a  complex  quantity  if  the  internal  fluid  acoustic  wavespeed  is  complex;  i.e., 
contains  a  loss  term  as  an  imaginary  quantity.  If  the  loss  factor  is  zero,  then  the  internal  fluid 
acoustic  wavespeed  is  purely  real,  and  y,  will  be  either  purely  real  or  purely  imaginary.  To  relate 
the  internal  acoustic  pressure  field  to  the  radial  shell  displacement  field,  conservation  of 
momentum  is  invoked  at  the  interface.  This  equation  is 


9 


(37) 


Pi 


d^u{x,a,t)  dpj{x,a,t) 


dP 


dr 


where  p,  is  the  density  of  the  interior  fluid  (kg  m'^).  Inserting  equations  (31)  and  (35)  into 
equation  (37)  allows  the  constant  A/ to  be  determined  and  the  pressure  field  to  be  written  as 


Pi{r)  = 


-CO  Pi  Joinr) 

Yi  Jl(//«) 


[Gi Ji  (/la)  +  G2 Yi  {Yxa)  +  G3 J, {y2a)  +  G4Y, (/2a)]  . 


(38) 


Inserting  equations  (31),  (32),  and  (38)  into  equation  (33)  yields  the  first  algebraic  boundary 
value  equation,  written  as 


a  Yi  Ji(/,«) 


+ 


+ 


(^11/2  ^^^13^2)^o(/2^) 


^12  ^ 

a  Yi  J|  (//«). 

^I2~^ll  I  Pi  JqC/i^) 

.  a  Yi  Ji  (//«), 

Cn  -g|i  ,  -(o' Pi  h^YP) 
a  Yi  Ji  (/,«). 


Y,(/,a) 


hiYiO) 


Y,(y,a) 


G,  =0. 


(39) 


The  second  boundary  condition  is  the  radial-longitudinal  shear  stress  in  the  shell  at  the  interface 
where  r  =  a  is  zero,  and  this  equation  is  written  as 


,  ,  ( du{x,a,t)  dw{x,a,t)\  . 

(7,^{x,a,t)  =  C44  — -■  + — ^  1  =  0  . 


\  dx 


dr 


(40) 


Inserting  equations  (31)  and  (32)  into  equation  (40)  yields  the  second  algebraic  boundary  value 
equation,  written  as 
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[c44(i^-r,^,)J,(r,a)]G, 

+  [c^(iA:-/,^,)Y,(/,o)]G, 
+[c^(\k-Y^4^)]^{y^a)]G^ 

+  [c^(^k-y,^,)Y,{y,a)]G,=^. 


(41) 


The  third  boundary  condition  is  a  force  balance  between  the  pressure  in  the  exterior  fluid,  an 
applied  radial  load,  and  the  normal  radial  stress  in  the  shell  at  the  interface  where  r  =  b.  This 
equation  is  written  as 


CJrr{x,b,t)^C]] 


du{x,b,t) 

dr 


+  C]2 


u{x,b,t)  dw{x,b,i) 


■  +  C]3 


dx 


=  Po{x,b,t)-  pg{x,t)  , 


(42) 


where  Pg(x,t)  is  an  applied  external  forcing  function  (N  m'^)  in  the  radial  direction  that  is 
assumed  to  be  at  a  discrete  wavenumber  and  frequency;  thus, 

Pe{x,t)  =  Pg  exp(ib:)exp(-i<y?)  ,  (43) 

and  p^(x,b,t)  is  the  scattered  acoustic  pressure  of  the  exterior  fluid  (N  m'^)  at  r=  b,  which 
satisfies  the  wave  equation  in  cylindrical  coordinates;  i.e.. 


d^Po(x,r,t)  ^  1  dpo{x,r,t)  ^  d^Po{x,r,t) _ 1  d^p„{x,r,t) 

dr^  r  dr  Sx^  ^  dt^ 

where  c„  is  the  acoustic  (or  compressional)  wavespeed  of  the  exterior  fluid  (m  s"').  Using  the 
infinite  length  of  the  cylinder  in  the  x:-direction  and  the  constraint  that  the  pressure  field  has 
vanished  when  r  approaches  infinity,  the  temporal  harmonic  solution  to  equation  (44)  is 

P„(x,r,t)  =  P„(A-)exp(ib:)exp(-i<yO  =  iVH[,'’(/„r)exp(ib:)exp(-i<yO  ,  (45) 

where  Hq  ^  denotes  a  zero-order  Hankel  function  of  the  first  kind  and 
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1/2 


Yo  = 


-k‘ 


K^o  J 


=  (^2-A:2)1/2 


(46) 


The  external  fluid  wavenumber  ■/„  can  be  a  complex  quantity  if  the  external  fluid  acoustic 
wavespeed  is  complex,  i.e.,  contains  a  loss  term  as  an  imaginary  quantity.  If  the  loss  factor  is 
zero,  then  the  external  fluid  acoustic  wavespeed  is  purely  real,  and  /o  will  be  either  purely  real  or 
purely  imaginary.  To  relate  the  external  acoustic  pressure  field  to  the  radial  shell  displacement 
field,  conservation  of  momentum  is  invoked  at  the  interface.  This  equation  is 


Po 


d  u{x,b,t)  dpg{x,b,t) 


dr 


dr 


(47) 


where  po  is  the  density  of  the  exterior  fluid  (kg  m'^).  Inserting  equations  (3 1 )  and  (45)  into 
equation  (47)  allows  the  constant  N  to  be  determined  and  the  pressure  field  to  be  written  as 


2  (0 

Po  (0  =  J 1  inb)  +  <^2 Y,  (7,6)  +  G3 J ,  {y^b)  +  G4 Y,  {y^b)]  .  (48) 

Yo  W^;\yoa) 

Inserting  equations  (31),  (32),  and  (48)  into  equation  (42)  yields  the  third  algebraic  boundary 
value  equation,  written  as 


(c„7,  +iA:c,3^,)Jo(7,6)  + 


Yo  ^TiYob), 


(c,,/,  +iA:c,3^,)Yo(/,6)  +  | 

(c„7,+i^c, 3^2)30(736)  + 
(c,,72+iA:c,3^2)Yo(726)  + 


-co^P.K\Yob) 

,  b  Yo  H1'>(7o6). 

'g|2-g|i  H^‘^(7^6) 


Yo  ^?{Yob) 


c,2-c„  H<‘'(7<,6) 


Yo  H{'>(7,6) 


Y,(7,6) 

J.(726) 

Y.(726) 


(49) 


^^4  = 


-P 
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The  fourth  boundary  condition  is  the  radial-longitudinal  shear  stress  in  the  shell  at  the  interface 
where  r  =  b  and  is  equal  to  an  applied  longitudinal  load  that  is  assumed  to  be  at  a  discrete 
wavenumber  and  frequency;  thus, 


^rx  0  ^44 


du{x,b,t)  ^  dw(x,b,t) 


dx 


dr 


(50) 


where  fgix,t)  is  an  applied  external  forcing  function  (N  m'  )  in  the  longitudinal  direction  that  is 
assumed  to  be  at  a  discrete  wavenumber  and  frequency;  thus, 

/^(x,/)  =  f;exp(iAx)exp(-i6;/)  .  (51) 


Inserting  equations  (31),  (32),  and  (51)  into  equation  (50)  yields  the  fourth  algebraic  boundary 
value  equation,  written  as 


[c^(i^-r,#,)j,(r,6)]G, 

+  [c^(iA:-r,^,)Y,(r,/j)]G, 

(52) 

-|-[C44(i^-/,^2)Jl(r2^)]^3 

+[c^0k-r,42mr2b)]G,=F,. 

Equations  (39),  (41),  (49),  and  (52)  are  now  written  in  matrix  form  as 

Ag  =  f,  (53) 

where  A  is  a  loiown  four-by-four  coefficient  matrix,  g  is  a  four-by-one  vector  that  contains  the 
four  unknown  wave  propagation  coefficients,  and  f  is  a  four-by-one  load  vector  that  represents 
the  external  forces  exciting  the  system.  (The  entries  of  the  matrix  and  vectors  in  equation  (53) 
are  given  in  appendix  A.)  The  wave  propagation  coefficients  are  now  found  by 


g  =  A  . 


(54) 


Once  the  wave  propagation  coefficients  are  known,  the  shell  displacements  can  be  calculated 
using  equations  (31)  and  (32),  the  exterior  pressure  field  can  be  calculated  using  equation  (48), 
and  the  interior  pressure  field  can  be  calculated  using  equation  (38). 
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3.  MODEL  VALIDATION 


The  model  is  now  validated  by  comparison  to  previously  developed  shell  theories.  First, 
the  fully  elastic  thick  shell  model  derived  in  section  2  is  compared  to  a  transversely  isotropic  thin 
shell  model.  From  a  previously  developed  isotropic  thin  shell  model,"  fluid  loading  is  added  to 
produce  a  longitudinal  equation  of  motion,  written  as 

. ho„E,  ^55^ 

dt^  Q-^rx^xr)  dx^  a{\-UrxOxr)  Sx 

and  a  radial  equation  of  motion,  written  as 


d^u(x,t) 
~dt 


=  -B 


d*u{x,t) 

dx* 


hE, 


u{x,t) 


hu^E^  dw{x,t) 


+  Pi  (a,x,  t)  -  (a,  X,  t)  -p^{x,t). 


where  B  is  the  flexural  stiffness  (N  m)  of  the  shell  and  is  given  by 


(56) 


B  = 


(57) 


Making  the  assumption  of  harmonic  response  in  space  and  time,  the  displacements  can  be 
written  as 

u(x,t)  =  U  exp(ifcc)exp(i<yr) ,  (58) 

and 


w(x,0  =  fFexp(ifcc)exp(i&»r)  .  (59) 

This  produces  the  matrix  equation 

Bu  =  p  ,  (60) 

where  the  unknown  displacements  are  contained  in  the  vector  u.  (The  entries  of  the  matrix  and 
vectors  in  equation  (60)  are  listed  in  appendix  A.)  The  unknown  displacements  are  determined 
using 


15 


(61) 


u  =  B  'p  . 

Once  the  displacements  are  known,  the  interior  pressure  field  for  this  model  can  be  calculated 
using 


Pi{r)  = 


-oP- Pi 

Yi  Jl  (//■«) 


and  the  exterior  pressure  filed  can  be  calculated  using 


(62) 


Poir)  = 


-oPp, 

Yo  H\^\roa) 


(63) 


Figure  2  is  a  plot  of  the  transfer  function  of  internal  pressure  at  r  =  0  divided  by  external 
forcing  function  in  the  radial  direction  versus  wavenumber.  Figure  3  is  a  plot  of  the  transfer 
function  of  internal  pressure  at  r  =  0  divided  by  external  forcing  function  in  the  longitudinal 
direction  versus  wavenumber.  In  both  figures,  the  upper  plot  is  the  magnitude  of  the  power 
expressed  in  the  decibel  scale  and  the  lower  plot  is  the  phase  angle  expressed  in  degrees.  The 
solid  line  is  the  transversely  isotropic  thick  shell  model  developed  in  section  2  and  the  dots 
correspond  to  the  transversely  isotropic  thin  shell  model  listed  as  equations  (55)  through  (63). 

In  this  example,  the  thickness  of  the  shell  was  small  (5.08  x  10"^  m)  and  the  frequency  was  low 
(100  Hz),  so  the  assumptions  of  the  thin  shell  model  are  valid  and  the  outputs  of  the  two  models 
should  reasonably  agree.  Figures  2  and  3  were  generated  with  the  following  parameters:  shell 
density  p=  1200  kg  m’^,  radial  Poisson’s  ratio  due  to  longitudinal  load  Vxr  =  0.48 
(dimensionless),  longitudinal  Young’s  modulus  Ex  =  2x  10^  N  m’^,  radial  Young’s  modulus 
Er  =  3  X  10*  N  m'^,  longitudinal  Poisson’s  ratio  due  to  radial  load  Orx  =  0.072  (dimensionless), 
shear  modulus  Gxr  =  6.76  x  10  N  m’  ,  inner  shell  radius  a  =  0.0759  m,  outer  shell  radius 
b  =  0.0765  m,  inner  fluid  density  =  800  kg  m’^,  inner  fluid  compressional  wavespeed 
c,  =  1300  m  s’',  outer  fluid  density  po  =  1000  kg  m'^,  and  outer  fluid  compressional 
wavespeed  Co  =  1500  m  s’'.  Based  on  these  values,  the  computed  stiffness  constants  are 
cii  =  3.15  X  10*  N  m’'',  ca  =  3.47  x  lO’  N  m’'',  cn  =  1.68  x  10*  N  m’'',  C33  =  2.16  x  10^  N  m’^ 
and  C44  =  6.76  x  10  N  m'  .  For  the  model  validation  problems  presented  here,  the  shell  has  zero 
damping;  however,  most  structures  have  some  loss  mechanism  associated  with  their  behavior. 
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Note  that  in  figures  2  and  3,  there  is  broad-based  agreement  between  the  thin  shell  model  and  the 
thick  shell  model.  It  is  noted  that  these  two  transfer  functions  are  of  interest,  and  these  specific 
outputs  will  be  investigated  in  the  remainder  of  this  report. 
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Figure  2.  Transfer  Function  of  Internal  Pressure  Divided  by  External  Radial  Excitation: 
Transversely  Isotropic  Thick  Shell  ( - )  and  Transversely  Isotropic  Thin  Shell  (•) 
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Figure  3.  Transfer  Function  of  Internal  Pressure  Divided  by  External  Longitudinal  Excitation: 
Transversely  Isotropic  Thick  Shell  ( - )  and  Transversely  Isotropic  Thin  Shell  (•) 
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Next,  the  fully  elastic  thick  shell  model  derived  in  section  2  is  compared  to  a  fully  elastic 
isotropic  thick  shell  model.  A  previously  developed  thick  shell  model^’^  is  based  on  Navier’s 
equation  of  motion  in  an  isotropic  solid,  written  in  vector  form  as 


(A  +  /J)VV  •  u  +  fuV^M  =  p 


(64) 


where  X  and  p  are  Lame  constants  (N  m'^)  of  the  shell  and  the  vector  u  represents  the 
displacement  field  (m).  This  equation  is  solved  and  coupled  to  the  inner  and  outer  pressure  field 
and  the  result  is  the  matrix  equation 

Cd  =  f  ,  (65) 

where  C  is  a  known  four-by-four  coefficient  matrix,  d  is  a  four-by-one  vector  that  contains  four 
unknown  wave  propagation  coefficients,  and  f  is  a  four-by-one  load  vector  that  represents  the 
external  forces  exciting  the  system.  (The  entries  of  the  matrix  and  vectors  in  equation  (65)  are 
given  in  appendix  A.)  The  wave  propagation  coefficients  are  now  found  by 


d  =  C"*f  .  (66) 

Once  the  wave  propagation  coefficients  are  known,  the  shell  displacements  can  be  calculated 
using 


u{x,  r,t)  =  [-DiOr  J,  (or)  -D^aY^  (or)  -  D^ikJ^(/3r)  -  DpkY^  (/?r)]exp(ifcc)  exp(-i6;/) ,  (67) 

and 

vv(x,  r,  t)  -  [D,iA:  Jo  (or)  -l-  D^iA:  ¥„  {ar)  +  {fir)  -l-  G^fiY^  (^6r)]  exp(ifcc)  exp(-i  cj?) .  (68) 

Once  the  displacements  are  known,  the  interior  pressure  field  for  this  model  can  be  calculated 
using 


P>ir)  =  («r)  -  Y,  {or)  -  D.xki,  (fir)  -  DpkY,  {Pr)] ,  (69) 
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and  the  exterior  pressure  filed  can  be  calculated  using 


Po ir)  = - — («'■) - Y, («r) -  J. (y9r) - (y9r)]  .  (70) 

Figure  4  is  a  plot  of  the  transfer  function  of  internal  pressure  at  r  =  0  divided  by  external 
forcing  function  in  the  radial  direction  versus  wavenumber.  Figure  5  is  a  plot  of  the  transfer 
function  of  internal  pressure  at  r  =  0  divided  by  external  forcing  function  in  the  longitudinal 
direction  versus  wavenumber.  In  both  figures,  the  upper  plot  is  the  magnitude  of  the  power 
expressed  in  the  decibel  scale  and  the  lower  plot  is  the  phase  angle  expressed  in  degrees.  The 
solid  line  is  the  transversely  isotropic  thick  shell  model  developed  in  section  2,  and  the  dots 
correspond  to  the  isotropic  thick  shell  model  listed  as  equations  (64)  through  (70).  In  this 
example,  the  transversely  isotropic  model  was  run  with  isotropic  material  properties,  so  the 
output  of  the  transversely  isotropic  model  should  reasonably  agree  with  the  output  of  the 
isotropic  model.  Figures  4  and  5  were  generated  with  the  following  parameters:  frequency 
/=  800  Hz,  shell  density  p  =  1200  kg  m'^,  radial  Poisson’s  ratio  due  to  longitudinal  load  Vxr  = 
0.48  (dimensionless),  longitudinal  Young’s  modulus  A  =  3  x  10*  N  m'^,  radial  Young’s  modulus 
Er  =  2)\  1 0*  N  m'^,  longitudinal  Poisson’s  ratio  due  to  radial  load  Urx  =  0.48  (dimensionless), 
shear  modulus  Gxr  =  1.01  x  10*  N  m'^,  inner  shell  radius  a  =  0.0762  m,  outer  shell  radius  h  - 
0.152  m,  inner  fluid  density  p-,  =  800  kg  m'^,  inner  fluid  compressional  wavespeed  c,  =  1300  m 
s'',  outer  fluid  density  po  =  1000  kg  m'^,  and  outer  fluid  compressional  wavespeed  Co  =  1500  m 
s  '.  Based  on  these  values,  the  computed  stiffness  constants  are  cn  =  A  +  2/y  =  2.64  x  10^  N  m'^, 
ci2  =  A  =  2.43  X  10^  N  m'^  C|3  =  A  =  2.43  x  10^  N  m'^  C33  =  A  +  2/2  =  2.64  x  10^  N  m'^  and  C44  = 
//=  1.01  X  10*  N  m'^.  Note  that  in  figures  4  and  5,  there  is  broad-based  agreement  between  the 
thin  shell  model  and  the  thick  shell  model. 
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Figure  4.  Transfer  Function  of  Internal  Pressure  Divided  by  External  Radial  Excitation: 
Transversely  Isotropic  Thick  Shell  ( - )  and  Isotropic  Thick  Shell  (•) 
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Figure  5.  Transfer  Function  of  Internal  Pressure  Divided  by  External  Longitudinal 
Excitation:  Transversely  Isotropic  Thick  Shell  ( - )  and  Isotropic  Thick  Shell  (•) 
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4.  HIGH  WAVENUMBER  APPROXIMATION 


Numerical  simulations  of  this  model  reveal  that  at  high  wavenumbers  the  A  matrix 
becomes  ill-conditioned  and  algorithmically  singular.  To  avoid  this  problem  in  any  analysis,  the 
outputs  of  the  model,  i.e.,  the  transfer  functions,  are  analyzed  differently  in  two  distinct 
regions — namely,  where  |A:|  <  5  /  a  and  where  |A:|  >  5  /  a .  In  the  region  where  |A:|  >  5  /  a ,  the 

model  outputs  are  calculated  in  such  a  manner  that  they  are  continuous  from  |A:|  <  5  /  a  to 
|A:|  >  5/a  and  they  are  proportional  to  Mk^.  This  is  written  in  equation  form  as 


Piir) 


Piir) 


and 


PXr) 


P^r) 


^  J  ^  e 


<1 

a 


W>1 


where 


Aq  = 


5 

aj 


Piir) 


and 


=  - 


5 

aj 


Piir) 


(71) 


(72) 


(73) 


(74) 


This  approximation  ensures  that  the  energy  falloff  past  the  flexural  wave  will  occur  in 
wavenumber  at  a  rate  that  is  observed  in  most  shell  models. 
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5.  NUMERICAL  EXAMPLE 


A  numerical  example  is  now  investigated  using  the  transversely  isotropic  shell  model 
derived  in  section  2.  Figure  6  is  an  image  of  the  magnitude  of  the  transfer  function  of  internal 
pressure  at  r  =  0  divided  by  external  forcing  function  in  the  radial  direction  versus  frequency  and 
wavenumber.  This  image  is  the  magnitude  of  the  power  expressed  in  the  decibel  scale  with  the 
scale’s  range  shown  as  a  colorbar  above  the  plot.  Figure  7  shows  constant  frequency  cuts  of 
figure  6  at  frequencies  of  500,  1000,  1500,  and  2000  Hz.  Figure  8  is  an  image  of  the  magnitude 
of  the  transfer  function  of  internal  pressure  at  r  =  0  divided  by  external  forcing  function  in  the 
longitudinal  direction  versus  frequency  and  wavenumber.  This  image  is  also  the  magnitude  of  the 
power  expressed  in  the  decibel  scale  with  the  scale’s  range  shown  as  a  colorbar  above  the  plot. 
Figure  9  shows  constant  frequency  cuts  of  figure  8  at  frequencies  of  500,  1000,  1500,  and 
2000  Hz.  Figures  6  through  9  were  generated  with  the  following  parameters:  shell  density  p  = 
1200  kg  m'^,  radial  Poisson’s  ratio  due  to  longitudinal  load  Vxr  -  0.48  (dimensionless), 
longitudinal  Young’s  modulus  Ex  =  2x  10\l-0.05i)  N  m'^,  radial  Young’s  modulus  Er  =  3  x 
10*(l-0.10i)  N  m longitudinal  Poisson’s  ratio  due  to  radial  load  Urx  -  0.0722(l-0.0498i) 
(dimensionless),  shear  modulus  Gxr  =  6.76  x  10\l-0.05i)  N  m'^,  inner  shell  radius  a  =  0.0762  m, 
outer  shell  radius  6  =  0.1 524  m,  inner  fluid  density  pi  =  800  kg  m'^,  inner  fluid  compressional 
wavespeed  c,  =  1300  m  s"',  outer  fluid  density  po  =  1000  kg  m'^,  and  outer  fluid  compressional 
wavespeed  Co  =  1 500  m  s"' .  Based  on  these  values,  the  computed  stiffness  constants  are  ci  i  =  3. 1 5 
X  10*(l-0.103i)Nm'^ci2  =  3.46x  10^(l-0.156i)  N  m■^  cu  =  1.68  x  10\l-0.108i)  N  m'^ 

C33  =  2.16  X  10^1-0. 0543i)  N  m'^,  and  C44  =  6.76  x  10\l-0.05i)  N  m'^.  For  this  problem,  the 
two  validation  models  used  in  section  3  are  not  capable  of  modeling  this  configuration;  thus,  no 
comparison  can  be  made  with  previously  available  solutions.  The  MATLAB  code  used  to 
generate  this  (and  the  previous)  example  is  included  as  appendix  B. 
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Figure  6.  Transfer  Function  of  Internal  Pressure  Divided  by  External  Radial  Excitation 
Versus  Wavenumber  and  Frequency  for  Numerical  Example 
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Figure  7.  Transfer  Function  of  Internal  Pressure  Divided  by  External  Radial  Excitation 
Versus  Wavenumber  for  (a)  500  Hz,  (b)  1000  Hz,  (c)  1500  Hz,  and  (d)  2000  Hz 
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Figure  8.  Transfer  Function  of  Internal  Pressure  Divided  by  External  Longitudinal 
Excitation  Versus  Wavenumber  and  Frequency  for  Numerical  Example 


26 


Magnitude  (dB)  Magnitude  (dB)  Magnitude  (dB)  Magnitude 


20 

0 

-20 

40 


/ 


\  „ 


-V 


\ 


(d) 


-X. 


-40  -30  -20  -10  0  10 

Wavenumber  (rad/m) 


20  30  40 


Figure  9.  Transfer  Function  of  Internal  Pressure  Divided  by  External  Longitudinal 
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6.  SUMMARY 


A  model  of  a  transversely  isotropic  thick  shell  with  fluid  loading  on  the  inner  and  outer 
surfaces  has  been  derived.  This  model  is  compared  to  two  previously  available  models  and  is 
shown  to  be  in  agreement  for  the  case  where  the  shell  is  transversely  isotropic  and  extremely  thin 
and  the  case  where  the  shell  is  isotropic  and  thick.  A  numerical  example  is  given  where  the  shell 
is  transversely  isotropic  and  thick.  A  calculation  to  bypass  the  high  wavenumber  instability  that 
is  typical  of  this  class  of  problems  is  included.  The  MATLAB  code  used  to  generate  the 
numerical  examples  is  also  included. 
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APPENDIX  A 

COEFFICIENTS  OF  MATRICES  AND  VECTORS 


This  appendix  contains  the  coefficients  of  the  matrices  and  vectors  from  the  models 
developed  in  this  report. 


The  entries  of  the  A  matrix  from  equation  (53)  are 
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The  g  vector  from  equation  (53)  is 

g  =  [Gi  G2  G3  ^4]^  .  (A-17) 


The  f  vector  from  equation  (53)  is 
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The  entries  of  the  B  matrix  from  equation  (60)  are 
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The  u  vector  from  equation  (60)  is 
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The  p  vector  from  equation  (60)  is 
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The  entries  of  the  C  matrix  from  equation  (65)  are 
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In  equations  (A-25)  through  (A-40),  the  constants  are  as  follows: 
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The  relationship  between  the  Lame  constants  and  the  material  properties  is 
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In  equations  (A-45)  and  (A-46),  Young’s  modulus  and  Poisson’s  ratio  can  be  along  any  axis 
because  the  material  is  isotropic. 
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APPENDIX  B 

MATLAB  SUBROUTINE  OF  MODEL 


%  - -Elast icShel ITF - Elastic  Shell  Transfer  Function 

% 

% This  program  produces  a  model  of  the  interior  pressure 

%  in  a  fluid  filled  shell  when  it  is  loaded  on  the  exterior 
%  by  a  normal  and  longitudinal  force.  The  shell  has  an 

%  outer  fluid  and  is  transversly  isotropic.  The  behavior  of 

%  the  shell  is  two-dimensional  fully  elastic.  Only  the  positve 

%  wavenumber  points  are  calculated  because  the  model  is  symmetric 

%  in  wavenumber . 

% 

%--Written  by  Andrew  J.  Hull  on  11/26/08 
% 

function  [  PiDPo,  PiDFo  ]  =  Elast icShellTF 
( f req, kmax, numpts , a, b, r , Ex, Er , nuxr , ro, roi , ci , roo, co) 

% 

% - Output  Variables 

%  PiDPo  =  Interior  Pressure  (at  r)  Divided  by  Exterior  Normal  Pressure 

%  PiDFo  =  Interior  Pressure  (at  r)  Divided  by  Exterior  Longituidnal  Force 

% 

% - Input  Variables 

%  freq  =  Frequency  (Hz) 

%  kmax  =  Maximum  wavenumber  (rad/m) 

%  numpts  =  Number  of  points  in  wavenumber 

%  a  =  Inner  shell  radius  (m) 

%  b  =  Outer  shell  radius  (m) 

%  r  =  Hydrophone  radius  (m) 

%  Ex  =  Modulus  in  the  axial  direction  (N/m^2) 

%  Er  =  Modulus  in  the  radial  direction  (N/m^2) 

%  nuxr  =  Poisson's  ratio  of  the  matrix  material  (dimensionless) 

%  ro  =  Density  of  the  shell  (kg/m^3) 

%  roi  =  Density  of  the  inner  fluid  (kg/m^3) 

%  ci  =  Wavespeed  of  the  inner  fluid  (m/s) 

%  roo  =  Density  of  the  outer  fluid  (kg/m^3) 

%  CO  =  Wavespeed  of  the  outer  fluid  (m/s) 

% 

% - - 

%-  Frequency  in  rad/s 
w  =  2  *  pi  *  freq; 

%" -Build  the  wavenumber  vector 

kvec  =  1  inspace  (  eps,  kmax,  numpts  )  ,* 

%--kahigh  is  the  wavenumber  cutoff  where  the  high  wavenumber 
%  approximation  is  used  in  the  analysis.  (This  can  be  changed) 
kahigh  =  5.0; 

%- -Determine  if  the  wavenumber  vector  has  to  be  broken  into  low 

%  and  high  regions 

if  (  a*kvec(end)  >  kahigh  ) 

start indexhigh  =  min  (  find  (  a*kvec  >  kahigh  )  ) ; 

klow  =  kvec ( 1 : 1 : start indexhigh-1 ) ; 
khigh  =  kvec (start indexhigh: 1 :end) ; 
else 

klow  =  kvec; 

end 

% 

%- -Transversely  isotropic  material  constants  from  physical  constants 
Gxr  =  Ex  /  {  2  *  (  1  +  nuxr  )  ) ; 

nurx  =  nuxr  *  (  Er  /  Ex  )  ; 

cll  =  Er  ♦  (l-nurx*nuxr)  /  (  (1+nurx)  *  (l-nurx-2*nurx*nuxr)  ) ; 
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cl2  =  Er  *  nurx  *  (1+nuxr)  /  (  (1+nurx)  *  (l-nurx-2*nurx*nuxr)  ) ; 

cl3  =  Er  *  nuxr  /  (l-nurx-2*nurx*nuxr) ; 

c33  =  Er  *  nuxr  *  (1-nurx)  /  (  nurx  *  (l-nurx-2*nurx*nuxr)  ) ; 

c44  =  Gxr; 

% 

% -  -  - 


%--Low  wavenumber  region  ka  <  kahigh 
b2  =  cll*c44; 

bl  =  (cll*c33  -  cl3^2  -  2*cl3*c44)  *klow.  ^2  -  (c44+cll)  ^'^ro*w^2  ; 

bO  =  ro'^2*w^4  -  ro*w^2*  (C33+C44 )  *klow.  ^2  +  c33*c44*klow.  ^4  ; 

gammal  =  sqrt  (  (  -bl  +  sqrt  (  bl.'^2  -  4*b0*b2  )  )  /  (  2  *  b2  )  )  ; 

gamma2  =  sqrt  (  (  -bl  -  sqrt  (  bl.^2  -  4*b0*b2  )  )  /  (  2  *  b2  )  ) ; 

zetal  =  (  ro*w^2  -  cll*gammal . ^2  -  c44*klow.^2  )  ./  (  i*klow . *gammal* (  cl3  + 

C44  )  )  ; 

zeta2  =  (  ro*w'^2  -  cll*gamma2 .  ^2  -  c44*klow.^2  )  ./  (  i*klow .  *gamma2*  (  cl3  + 

c44  )  )  ; 

% 

%-- Inner  fluid  load 

gammai  =  sqrt  (  (w/ci)  "^2  -  klow.'^2  ); 
gamma i  =  gammai  +  (  gammai  ==  0  )*eps; 

fluidiload  =  (  (-w'^2*roi)  ./  gammai  )  .*  (  bessel j  ( 0 , gammai*a)  ./ 

besselj (l,gammai*a)  ) ; 


%  -Srr(a)  =  inner  fluid  load; 

Amatll  =  (  cll  *  gammai  +  i  *  klow  *  cl3  zetal  ) 


(  1  /  a  ) 


fluidiload 


(  (  cl2  -  cll  ) 

besselj (l,gammal*a)  ; 

Amatl2  =  (  cll  *  gammai  +  i  *  klow  *  cl3  .*  zetal  ) 


(  1  /  a  ) 


fluidiload 


(  (  cl2  -  cll  ) 

bessely ( 1 , gammai *a)  ; 

Amatl3  =  (  cll  *  gamma2  +  i  *  klow  *  cl3  zeta2  ) 


(  1  /  a  ) 


fluidiload 


(  (  cl2  -  cll  ) 

besselj (l,gamma2*a)  ; 

Amatl4  =  (  cll  *  gamma2  +  i  *  klow  *  cl3  .*  zeta2  ) 


(  (  cl2  -  cll  ) 

bessely (1 , gamma2*a) ; 

% 

%--Srz(a)  =  0; 

Amat21  =  c44  * 


(  1  /  a  ) 


fluidiload 


*  besselj (0,gammal*a) 

★ 

*  bessely (0 , gammal*a) 

* 

*  besselj (0,gamma2*a) 

* 

*  bessely (0 , gamma2*a) 


i 

*  klow 

-  gammai 

,*  zetal  ) 

.*  besselj ( 1 , gammal*a) 

i 

*  klow 

-  gammai 

.  *  zetal  ) 

. *  bessely ( 1 , gammal*a) 

i 

*  klow 

-  gamma 2 

.*  zeta2  ) 

.*  besselj ( 1 , gamma2*a) 

i 

*  klow 

-  gamma 2 

.  *  zeta2  ) 

.*  bessely ( 1 , gamma2*a) 

%- -Outer  fluid  load 

gammao  =  sqrt  (  (w/co)^2  -  klow. ^2  ); 
gammao  =  gammao  +  (  gammao  ==  0  )*eps; 
fluidoload  =  (  (-w'^2*roo)  ./  gammao  )  .* 

besselh ( 1 , 1 , gammao*b)  ); 

% 

%--Srr(b)  =  outer  fluid  load 
Amat31  =  (cll  *  gammai  +  i  *  klow  *  cl 3 
+  ... 

(  (  cl2  -  cll  )  *  (  1  /  b  ) 

besselj (l,gammal*b) ; 

Amat32  =  (  cll  *  gammai  +  i  *  klow  *  cl3 
+  ... 

(  (  cl2  -  cll  )  *  (  1  /  b  ) 

bessely ( 1 , gammal*b) ; 

Amat33  =  (  cll  *  gamma2  +  i  *  klow  *  cl3 
+  ... 


(  besselh (0,1, gammao*b)  . / 


,  *  zetal  ) 
fluidoload 
,  *  zetal  ) 
fluidoload 
,  *  zeta2  ) 


.*  besselj ( 0 , gammal*b) 


.  *  bessely (0, gamma l*b) 


.*  besselj (0,gamma2*b) 
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(  (  cl2  -  cll  )  *  (  1  /  b  )  -  fluidoload  )  .* 

besselj { 1 , gamma2*b) ; 

Amat34  =  (  cll  *  gamma2  +  i  *  klow  *  cl3  .*  zeta2  )  .*  bessely ( 0 , gamma2*b) 

+  ... 

{  {  cl2  -  cll  )  *  (  1  /  b  )  -  fluidoload  )  .* 

bessely (1 , gamma2*b) ; 


%--srz (b) 
Amat 41  = 

C44 

* 

( 

i 

* 

klow 

-  gamma 1 

Amat42  = 

C44 

* 

( 

i 

* 

klow 

-  gamma 1 

Amat43  = 

C44 

* 

( 

i 

* 

klow 

-  gamma2 

Amat 4 4  = 

% 

C44 

* 

( 

i 

* 

klow 

-  gamma 2 

zetal  )  .*  besselj ( 1 , gammal*b) ; 
zetal  )  .*  bessely ( 1 , gammal*b) / 
zeta2  )  .*  besselj ( 1 , gamma2*b) ; 
zeta2  )  .*  bessely (1 , gamma2*b) ; 


%- -Determinant  of  Amat 
DetA  =  Amatll . *Amat22 . *Amat33 . *Amat44 
-Amatll . *Amat32 . *Amat23 . *Amat44 
Amatll . * Amat 4 2 . * Amat 2 3 . * Amat 3 4 
-Amat21 . *Amatl2 . *Amat33 . *Amat44 
Amat21 . *Amat32 . *Amatl3 . *Amat44 
-Amat21 . *Amat42 . *Amatl3 . *Amat34 
Amat 31 . *Amatl2 . * Amat 2 3 . * Amat 4 4 
-Amat31 . *Amat22 . *Amatl3 . *Amat44 
Amat 31 . * Amat 4 2 . *Amatl3 . *Amat24 
-Amat 41 . *Amatl2 . * Amat 2 3 . *Amat34 
Amat41 . *Amat22 . *Amatl3 . *Amat34 
-Amat 41 . *Amat32 . *Amatl3 . * Amat 2 4 


-  Amatll . * Amat 2 2 . *Amat34 . * Amat 4 3  + 
+  Amatll , *Amat32 . *Amat24 . *Amat43  + 

-  Amatll . *Amat42 . *Amat24 . *Amat33  + 
+  Amat21 . *Amatl2 . *Amat34 . *Amat43  + 

-  Amat21 . *Amat32 . *Amatl4 . *Amat43  + 
+  Amat21 . *Amat42 . *Amatl4 . *Amat33  + 

-  Amat31 . *Amatl2 . *Amat24 . *Amat43  + 
+  Amat31 , *Amat22 , *Amatl4 . *Amat43  + 

-  Amat 31 , *Amat42 , * Amat 14 . * Amat 2 3  + 
+  Amat41 , *Amatl2 . *Amat24 . *Amat33  + 

-  Amat41 , *Amat22 . *Amatl4 . *Amat33  + 
+  Amat4 1 , *Amat32 . *Amatl4 , *Amat23 ; 


%-  Protect  against  a  zero  divide 
DetA  =  DetA  +  (  DetA  ==  0  )*eps; 

% 

%- -Inverse  terms  of  A  (Not  all  terms 
invA13  =  (  Amatl2 . *Amat23 . *Amat44  - 
-Amat22 . *Amatl3 . *Amat44  + 
Amat42 . *Amatl3 . *Amat24  - 

% 

invA14  =  (  -Amatl2 . *Amat23 . *Amat34  + 
Amat22 . *Amatl3 . *Amat34  - 
-Amat32 . *Amatl3 . *Amat24  + 

% 

% 

invA23  =  (  -Amat23 . *Amat44 . *Amatll  + 
-Amat41 . *Amatl3 . *Amat24  + 
Amat 21 . *Amatl3 . * Amat 4 4  - 

% 

invA24  =  (  Amat23 . *Amat34 . *Amatll  - 
Amat31 . *Amatl3 . *Amat24  - 
-Amat21 . *Amatl3 . *Amat34  + 


are  needed) 

Amat 12 . * Amat 2 4 . * Amat 4 3 
Amat22 . *Amatl4 . *Amat43 
Amat42 . *Amatl4 . *Amat23 

Amat 12 . *Amat24 . *Amat33 
Amat22 . *Amatl4 . *Amat33 
Amat32 . *Amatl4 , *Amat23 


Amat24 . *Amat43 . *Amatll 
Amat41 . *Amatl4 . *Amat23 
Amat21 . *Amatl4 . *Amat43 

Amat24 . *Amat33 . *Amatll 
Amat 31 . *Amatl4 . *Amat23 
Amat21 . *Amatl4 . *Amat33 


+ 

+ 

+ 

+ 


+ 

+ 


+ 

+ 


./  DetA; 

./  DetA; 


./  DetA; 

./  DetA; 


invA3  3 


invA34 


% 

invA4  3 


invA4  4 


(  Amat22 . *Amat44 . *Amatll 
Amat41 . *Amatl2 . *Amat24 
-Amat 21 . * Amat 12 . *Amat44 


-  Amat24 . *Amat42 . *Amatll  + 

-  Amat4 1 . *Amatl4 . *Amat22  + 
+  Amat21 . *Amatl4 . *Amat42  ) 


'j 


DetA; 


(  -Amat22 . *Amat34 . *Amatll 
-Amat 31 . * Amat 12 . *Amat24 
Amat 21 . * Amat 12 . * Amat 3 4 


+  Amat24 . *Amat32 . *Amatll  + 
+  Amat31 . *Amatl4 . *Amat22  + 

-  Amat21 . * Amat 14 . * Amat 3 2  ) 


DetA; 


(  -Amat22 . *Amat43 . *Amatll 
-Amat41 . *Amatl2 . *Amat23 
Amat 21 . * Amat 12 . * Amat 4 3 


+  Amat23 . *Amat42 . *Amatll  + 
+  Amat4 1 . *Amatl3 . *Amat22  + 

-  Amat21 . *Amatl3 . *Amat42  ) 


./  DetA; 


(  Amat22 . *Amat33 . *Amatll 
Amat 31 . *Amatl2 . * Amat 2 3 


Amat23 . *Amat32 . *Amatll  +  ... 
Amat31 . *Amatl3 . *Amat22  +  . . . 


B-3 


-Amat21 .  *Amatl2  .  *Amat33  +  Ainat21 .  *Amatl3  .  *Anat32  )  ./  DetA; 

% 

%- -Shell  displacement  radial  direction  at  a  due  to  ext6?rnal  radial  pressure 
ShellDispRadDPo  =  -invA13 . *bessel j (1 , gammal*a)  -  invA33 . *bessel j ( 1 , gamma2*a)  + 

-invA23 . *bessely (1 , gammal*a)  -  invA43 . *bessely ( 1 , gamma2*a) ; 

% 

%-- Interior  fluid  pressure  at  r  due  to  external  radial  pressure 

PiDPo  =  (  (-w^2*roi)  ./  gammai  )  .*  ShellDispRadDPo  .*  (  bessel j ( 0 , gammai*r ) 

./  besselj (l,gammai*a)  ); 

% 

%- -Shell  displacement  radial  direction  at  a  due  to  external  longituidnal  force 
ShellDispRadDPo  =  invA14 . *bessel j ( 1 , gammal*a)  +  invA34 . *bessel j ( 1 , gamma2*a)  + 

invA24 . *bessely (1, gammal*a)  +  invA44 . *bessely ( 1 , gamma2*a) ; 

% 

%-- Interior  fluid  pressure  at  r  due  to  external  longitudinal  force 

PiDFo  =  (  (-w^2*roi)  ./  gammai  )  .*  ShellDispRadDPo  .*  (  bessel j ( 0 , gammai*r) 

./  bessel j (1 , gammai*a)  ); 

% 

% - 

%--High  wavenumber  region  ka  >=  kahigh 
if  (  exist (' khigh ' )  ==  1  ) 

Azero  =  PiDPo (end)  *  (klow (end) ) ^4 ; 

Bzero  =  PiDPo (end)  *  (klow (end) ) ^4 ; 

% 

PiDPohigh  =  Azero  ./  (  khigh. ^4  ); 

PiDPohigh  =  Bzero  ./  (  khigh. ^4  ); 

% 

PiDPo (  max ( size (klow) +1)  :  max(size (kvec) )  )  =  PiDPohigh; 

PiDPo (  max (size (klow) +1)  :  max ( size (kvec) )  )  =  PiDPohigh; 

% 

end 

% 

% - 

%-- Populate  the  output  vectors  with  the  negative  wavenumber  response 
PiDPo  =  [  f liplr (PiDPo (2 : end) )  PiDPo  ]; 

PiDPo  =  [  fliplr (PiDPo(2:end) )  PiDPo  ]; 

% 

end 

% -  - 
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